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In a recent article (Forterre & Pouliquen 2001), we have reported a new instability ob- 



served in rapid granular flows down inclined planes that leads to the spontaneous for- 
mation of longitudinal vortices. From the experimental observations, we have proposed 
an instability mechanism based on the coupling between the flow and the granular tem- 
perature in rapid granular flows. In order to investigate the relevance of the proposed 
mechanism, we perform in the present paper a three-dimensional linear stability anal- 
ysis of steady uniform flows down inclined planes using the kinetic theory of granular 
flows. We show that in a wide range of parameters, steady uniform flows are unstable 
under transverse perturbations. The structure of the unstable modes are in qualitative 
agreement with the experimental observations. This theoretical analysis shows that the 
kinetic theory is able to capture the formation of longitudinal vortices and validates the 
instability mechanism. 



1. Introduction 

Unlike classical fluids which are well-described by Navier-Stokes constitutive equa- 
tions, granular flows still lack an unifying description. For slow deformations at high 
density, multi-body interactions and friction between grains control the dynamics of 
the flow. On the other hand, when the energy injected into the material is large, the 
particles are strongly agitated and interact mainly via instantaneous collisions. In this 
coUisional regime, the material can be compared to a gas and a "gran ular temperature" 
can be defined in relation with the random motion of the particles ( |Ogawa et al 1980 ; 
Campbell 1990 ). This analogy between a coUisional granular flow and a molecular gas 
has led to the development of a kinetic theory for rap id granular flows (Jenki ns & Sav- 
age 1983; |Haff 1983[ |Lun et al 19"84[ |Brey et al 1998|; ^ela fc Goldhirsch 1998D which is 
inspired by the kinetic theory of dense gas ( Chapman fc Cowling 1970 ). However, the 
main difference with classical molecular gases is that the collisions between granular parti- 
cles are inelastic. If no energy is supplied to the system, the granular temperature decays 
rapidly because each collision removes kinetic energy from the particles. It can be shown 
that a free dissipative gas can form dense clusters and eventually collapses in a finite time, 
as a consequence of inelasticity ( McNamara fc Young 199^ ; Goldhirsch & Zanetti 1993). 
In order to maintain the coUisional regime, energy must therefore be supplied to the sys- 
tem, for instance by strongly shaking the boundaries (Warr et al 1995; Falcon et al 1999; 



Rouyer & Menon 2000). In a rapid granular flow, another way to maintain granular tem- 



perature is to impose a shear to the mean flow. In that case, the granular temperature 
results from a balance between the production by the shear work and the lost due to 
the inelastic collisions. In return, the granular temperature influences the mean flow 
because the pressure and the transport coefficients (e.g. viscosity, thermal conductivity) 
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depend on the temperature, as in a molecular gas. The coupling between the granular 
temperature and the mean flow is one of the fundamental properties of rapid granular 
flows. Understanding the role of this coupling is thus important in order to better define 
the analogies and the differences between classical fluid and granular flows. 

Recently, we have reported an instability observed in rapid granular flows down rough 
inclined planes which seems to result from this coupling (Forterre & Pouliquen 2001). 
The experimental set-up was a rough inclined plane as sketched in figure |l|(a). For high 
inclinations and large openings of the reservoir, the free surface deforms in a very regular 
pattern of longitudinal streaks parallel to the flow direction (see the picture in figure |l|a). 
Velocity measurements of the grains have revealed that the streaks correspond to the 
formation of longitudinal vortices as sketched in figur e ^(b). Alt hough such structures 
are common in fluid mec hanics (e.g. G rtler vortices, see Baric 1994 or streaks in turbulent 



boundary layers, see Kachanov 1994 ), they had not been observed before in granular 



flows. The main experimental observations are the following: 

• The wavelength A of the surface deformation scales with the mean thickness h of 
the flow (A ~ 2 - 3/i). 

• The troughs correspond to downward part of the flow while the crests correspond 
to upward part of the flow (sec figure |l|b) . 

• The a;-component of velocity is greater in the troughs than in the crests. 

• Troughs are dense whereas crest are dilute (see figure |^). 

• The pattern drifts slowly in the transverse direction. 

These observations suggest the following instability mechanism to explain the formation 
of the longitudinal vortices (Forterre & Pouliquen 2001). Because of the collisions with 
the rough bottom, particles close to the plane are strongly agitated i.e. the granular 
temperature is high at the bottom. Since high temperature means low density, density 
may become lower at the bottom than at the free surface. The flow is then mechanically 
unstable under gravity yielding longitudinal vortices. The mechanism we propose is 
analogous to the classical Rayleigh-Benard instability observed when a fluid is heated 
from below. However, in the present case the heating is not imposed by a thermostat but 
is created by the flow itself through the coupling between temperature and shear specific 
to granular media. 

The above explanation for the vortices formation is based on density profile inver- 
sion. In a granular dissipative gas, the density profile actually results from a complex 
balance between gravity, collisions and dissipation and its prediction is not straightfor- 
ward. Density profiles with higher density at the free surface than below have been 
observed in numerical simulations of rapid granular flows using discrete element meth- 
ods ( Campbell fc Brenncn 1985 ; Azanza 1998 ). However, no instability was observed 
because the simulations were two-dimensional. In the present paper, we want to inves- 
tigate the relevance of the proposed instability mechanism by using the kinetic theory 
of granular flows. We present a linear stability analysis of steady uniform flows down 
inclined planes, in the framework of the classical kinetic theory of Lun et al (1984). 
This theory provides a set of hydrodynamic equations coupling the density, the veloc- 
ity and the granular temperature under the assumption of instantaneous binary inelas- 
tic collisions. Although steady uniform flows down inclined planes have been studied 
within this framework ( Anderson fc Jackson 199"^ ; Ahn et al 1992; Azanza et al 1999), 
no stability analysis has been performed. Stability studies of rapid granular flows with 



the kinetic theory have mainly focused on two-dimensional Couette flows ( Ravage 1992 
Chi-Hwa Wang et al 199% Alam & Nott 1998; |Nott et al l"999| ). 

We do not expect from the present study a complete and quantitative description 
of the longitudinal vortices instability, since the applicability of the kinetic theory is 
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Figure 1. (a) Experimental set-up. The picture is a top view of the free surface lit from side 
when the instability is fully developed (the granular material is sand 0.25 mm in mean diameter, 
the angle of inclination 9 of the plane is 6 — 41°, the opening hg of the reservoir is hg = 13 
mm), (b) Sketch of the flow deduced from the measurements showing the surface deformation, 
the longitudinal vortices and the density variations. 



known to be limited. The lack of separation between the microscopic and macroscopic 
scales inherent to inelastic gases (Tan & Goldhirsch 1999) and the existence of multibody 
interactions when density is high, are serious difficulties which are the subject of active 
researches (Goldhirsch 1999). However, the kinetic theory qualitatively contains the 
most important character of rapid granular flow i.e. the coupling between shear flow and 
particle agitation due to the inelasticity of collisions. This coupling being the core of the 
instability mechanism we propose, this analysis should capture at least qualitatively the 
formation of longitudinal vortices. 

The paper is structured as follow. The Section 2 gives the governing equations and 
the boundary conditions we will use in the paper. Basic flows i.e. steady uniform flows 
down inclined planes are studied in Section 3. Section 4 is devoted to the linearization 
of the equations around the basic state and the numerical method. The results of the 
stability analysis are given in Section 5. Comparison with the experiment and discussion 
are presented in Section 6. Concluding remarks are given in Section 7. 



2. Kinetic theory of granular flows 

In this section, we recall the equations of the kinetic theory of granular flow we will 
use to investigate the formation of longitudinal vortices. 

2.1. Governing equations 

The kinetic theory of granular flows provides a set of hydrodynamic equations coupling 
the density p, the mean velocity u, and the granular temperature T under the assumption 
of instantaneous binary inelastic collisions (the granular temperature T is defined by 
^ < (5u^ >, where Su is the random velocity fluctuations). In the presence of gravity, 
the hydrodynamic equations are: 



(2.1) 
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P^=Pg + V.S, (2.2) 

= S:Vu-V.q-7, (2.3) 

where A/At — d / dt+u.V . The first equation ( pTl| ) is the conservation of mass. The second 
equation (p]2) comes from the conservation of momentum, where S is the stress tensor 



and g is the gravity acceleration. The third equation (2.3) is the energy equation where 
the temporal variation of the random kinetic energy is balanced by three terms. The 
term S : Vu represents the production of fluctuation energy due to the work of the stress 
S during shear. The term — V.q, where q is the flux of fluctuation energy, represents the 
conduction term. The term 7 is the collisional rate of energy dissipation. The distinctive 
feature of rapid granular flows lies in the dissipative term 7 due to the inelastic collisions 
between particles. This last term is responsible for the coupling between the mean flow 
and the granular temperature. 

The kinetic theory gives the constitutive relations for S, q and 7 as a function of the 
flow variables p, u and T. For the present purpose, we will use the closure due to Lun et 
al (1984). The total stress tensor S, the heat flux q and the rate of energy dissipation 7 
are written as: 

S = - {P{u, T) - T)V.u} I + 27^{v, T)S, (2.4) 
q= -X(z/,r)VT, (2.5) 

l^^{l~e')h{v)Ti. (2.6) 

Here I is the identity matrix, S ~ \{uij + Uji) — is the deviatoric part of the 

rate of deformation, pp is the particle density, v — p/ pp is the solid fraction and d is the 
particle diameter. We have omitted in ( |2.5| ) the contribution of the gradient of compacity 
to the heat flux. This term increases considerably the algebra while it has little effect on 
the qualitative nature of the results. Note that the collisional rate of energy dissipation 
7 is proportional to (1 — e^), where e is the coefficient of inelasticity of the particles 
(0 < e < 1). As in classical dense gases, the pressure P{v,T), the viscosities {rj{v,T), 
^(z^, T)) and the thermal conductivity K{v,T) depend on the solid fraction v and the 
temperature T: 

P{y,T) = PpfMT, 

n{v,T) = ppdf2{iy)TK , . 

a^,T) = Ppdfo{>^)Ti, ^ ■ ' 

K{iy,T) = ppdh{v)T^-, 

where the dimensionless functions fo{v)^ /i(^); hiv), f3,{y) and f^iv) are given in table 
|l| ( p!^un et al 1984 ). These functions involve the radial distribution function go{v)- In the 



following, we shall used the one suggested by Lun and Savage (1986): 

ffo(^) - \t^^ (2-8) 

1 - 



where Vm is the maximum solid fraction = 0.6 in the following). This radial dis- 
tribution function is suitable for free surface flows since the resulting equations have no 
singularity at 1/ = 0. We have checked that this choice does not change qualitatively the 
results by using other kind of radial distribution functions, such as the Carnahan-Stirling 



radial distribution (Jenkins & Savage 1983). Finally, all the equations that follow will 
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Table 1. Dimensionless constitutive functions. 77 = ^(1 + e) 



be written in term of non-dimensional variables defined by: 



{x,y,z) = -(a;,y,z), t 



■ t, v 



Pp 



1 u, f=^r. (2.9) 
gd gd 



For sake of simplicity in the notation, the tildes will be omitted and the solid fraction v 
will be called density in the following . 

2.2. Boundary conditions 

In order to solve the problem of granular flows down rough inclined planes, we have to 
specify boundary conditions for v, u and T both at the free surface of the flow and at 
the plane. Unlike classical fluid, the velocity in rapid granular flows does not vanish at 
a fixed solid boundary. Therefore, the rough plane may act as a source (resp. a sink) of 
fluctuating energy whether the shear work of the slip velocity is larger (resp. smaller) 
than the inelastic loss due to collisions with the plane. Boundary conditions for rapid 
granular flows at a rough surface have been the subject of extensive researches (Hui et 
al 1984; Jenkins & Richman 1986; Johnson et al 1990; Goldhirsch 1999). Here we will 
use the heuristic approach of Johnson et al (1990) relating the tangential stresses t.S.n 
and the heat flux q.n at the plane to density, slip velocity Ug and temperature by: 



t.S.n = rj*{iy,T) ||us|| , at the plane, 
q.n = (u.S).n — 7*(i^, r), at the plane. 



(2.10) 
(2.11) 



The vector t is parallel to the plane and defined as t = Ug/ ||us|| . The vector n is normal 
to the plane. Finally the function r]*{i', T) and 7*(j^, T) are given by: 



r7*(i.,r) = ^feHT^ 



(2.12) 



where the dimensionless functions f^iv) and //(j^) are given in table l[ The relation 
( ^.10 ) is a transfer of momentum balance at the plane. Equation (2.11) expresses that 
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heat can be produced at the plane if the shear work is stronger than the loss of energy 
due to collisions with the plane. The boundary conditions (2.1C) and (2.11) depend 
on two diniensionless parameters, 4> and e^,, which are related to the wall properties. 
The parameter is the particle-wall coefficient of restitution. In the following, 
will be taken equal to the particle-particle coefficient of restitution e. The parameter 
<f) is related to the rate of momentum transfer to the flow by the collision with the 
plane. Its value can be connected to the wall geometry in the case of 2D-flows of disk 
(Jenkins & Richman 1986). For a rough plane made of close-packed disks, one obtains 
a value </> ~ 0.1. For 3D-flows, one can expect lower values of </> since collisions do not 
always occur in the shear plane. We will use in the following the value cj) ~ 0.05 and will 
discuss its influence later. 

At the free surface, the stresses and the energy flux must vanish. However, the loca- 
tion of the free surface is not known a priori and its definition is somewhat artificial for 
very dilute flows. Rather than define the location of the free surface, we impose stresses 
and energy flux to vanish when the distance from the plane goes to infinity. It should 
be notice that the boundary conditions used here are somewhat different than those 



used in previous studies (Johnson et al 1990; Anderson & Jackson 1992; Ahn et al 1992 
In the work of Anderson 



Azanza et al 1999). 



granular layer is taken as a control parameter. 



Jackson (1992), the thickness h of the 
Numerical difficulties then arise when 

Ahn et 



trying to match the stress-free condition at the free surface (Johnson et al 199C) 
al (1992) do not define the free surface but impose stresses to vanish at infinity. However, 
they arbitrarily fixe the density, velocity and granular temperature at the plane. The 
same procedure is used by Azanza et al (1999). Here, we adopt a mixed point of view 
since at the plane the boundary conditions of Johnson et al (1990) are satisfied whereas 
at infinity the procedure of Ahn et al (1992) is chosen. 



3. Steady uniform flows 

The first step in order to perform a linear stability analysis is to determine the basic 
flow i.e. two-dimensional steady uniform flows. Steady uniform flows down inclined 
planes have already been studied in the framework of the kinetic theory (Anderson & 
Jackson 1992; Ahn et al 1992; Azanza et al 1999| ). It is not in the scope of the present 
study to make an extensive investigation of steady uniform flows. Rather, we shall focus 
on the shape of the density profile, which plays an important role in the instability 
mechanism we propose. 



3.1. Equations for steady uniform flows 

Sl) with boundary conditions ( 2.1C ) and ( 2.11 ) to two- 



We apply equations (|2.1| 
dimensional steady uniform flows down inclined planes 
density, velocity and temperature in the following form: 



We thus look for solution for 



iy{x,y,z,t) 
u(x,y,z,t) 
T{x,y,z,t) 



fo(z), 
uo(z)x, 
To{z). 



(3.1) 



In such a flow the derivatives parallel to the plane are zero and the mass-conservation 
equation (2.1) is satisfied. The momentum-conservation equation (^.2[) in the flow direc- 
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tion (x-dircction) and in the direction normal to the plane (z-direction) becomes: 

dS. 



dz 
dSz^o 



dz 



= —vq smt 



vq cos 6*, 



(3.2) 
(3.3) 



where 9 is the angle of inclination of the plane, Tixzo = /2(^'o)To^^^duo/dz and Szzq 
— /i(i'o)'7o- The energy-conservation (2_^) simplifies to: 







duQ dq^Q 



7o, 



(3.4) 



where = -/sMToi/^dTo/dz and 70 = (1 - e^)f5MTo^^^- 

From ( p.2| ) and (3^) together with boundary conditions ^ xzo — '^zzo_^^ ^it infinity, one 
obtains S^zo = — tanO Ti^zq- Therefore, equations (3.2), (3.3) and (3.4) can be written 
in the following form: 



dvQ 
dz 

dup 
dz 

d^To 
dz^ 



1^0 cos ( 



- T^tan^To^ 

/2(i^0) 



[l - e^)f5Mf2M - fiHi^o)tan0 



f2{vo)f3{j^o) 



Tq H : — : — r, — ^ — cos 0— — 

/3(^^o)/i'(^o)To dz 



1 



dTo 

dz 



(3.5) 
(3.6) 

(3.7) 



JaMfi {vq) 

where the dashes in fi and /s' denote differentiation with respect to vq. 
This system of three nonlinear ordinary differential equations requires boundary condi- 
tions both at the plane and at infinity. At infinity, the stresses and the flux of energy 
must vanish. It can be shown ( Ahn ct al 1992 ) that these conditions are satisfied if and 
only if the temperature gradient vanishes at infinity: 

dio 
dz 



when 



(3.8) 



At the plane, the relations (2.1C) and (2.11) become: 

0/6(i^o)7o^uo at z = 0, 

{l-e^)f,(v^T^ 



f , 1 duo 
dz 



-fsMTo' 



dTo 



dz 



f f 1 duo 
dz 



at z 



0. 



(3.9) 
(3.10) 



These relations can be combined with equation ( |3.6| ) to express the temperature and the 
gradient of temperature at the plane as a function of density and velocity at the plane: 



To 

dTo 
dz 



/i'(i.o)tan2 
h{vo) 



uo at z = 0, 

, 0/6(^o)(l-e^)/7(^o) 
/i2(,.„)tan20 



Uo 



at z = 0. 



(3.11) 
(3.12) 



In order to integrate numerically the above boundary value problem, a "shooting method" 
has been chosen. The angle 8 and the density at the wall i'o(O) are the two control 
parameters whereas the slip velocity Us — uo(0) is taken as the shooting parameter. Then, 
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Figure 2. Phase diagram for steady uniform flows in the (h, v) plane. The thin solid lines give 
the contours of constant angle 6. The grey zone is the domain of non-inverted density profiles, 
(e = 0.6, (j> = 0.05). 



t he te mper ature and the gradient of temperat ure a t th e plane are given by the equations 
( ^.11[ ) and ( |3.12| ). The differential equations (|U)-(|3!|) can therefore be integrated from 



z = to z = oo. This procedure is repeated for different values of the shooting parameter 
Us = ■tio(O) until the boundary condition at infinity dTo/dz = is satisfied. In practice, 
a fourth-order Runge-Kutta method is employed and the equations are integrated up 
to a maximum value z = Zmax- The boundary condition at infinity is considered to be 
satisfied when the gradient of temperature is smaller than 10^^ at z = Zmax- 

3.2. Results 

By applying the above procedure to different values of the angle 9 and the density i'o(O), 
we were able to explore the whole set of steady uniform flows. However, in order to 
help the physical interpretation of the results, we present the results in the {h, v) plane 
where h is the thickness of the flow and v the mean density. The thickness h is deflned a 
posteriori by the value z = h for which the density is 1% of the maximum density inside 
the flow. The mean density i> is simply given hy v = (l/h) vo{z)dz. The relation 
between {6, i^o(O)) and {h, v) is univocal i.e. a single solution exists for a given set 
(ft,, v) (which is not the case if one chooses the flow rate as a control parameter, see 
Johnson et al 199C| and Anderson & Jackson 1992). The relation between (0, t'o(O)) and 



[h, v) is shown in figure |^ for e = 0.6 and (p = 0.05. In this figure, the thin solid lines 
are obtained for constant inclination 9 by varying vq{Q). We observe that the angle 
has a strong influence on the thickness of the flow h, the thickness of the flow rapidly 
increases when increasing the angle. For angles 9 > 23°, the thickness diverges whereas 
for 9 < 14°, the thickness becomes less than one grain diameter. Therefore, for a given 
set of parameters (e, (j>) there exists a flnite range of angles where steady uniform flows 
are obtained. These results are consi stent with analytical solut ions for steady uniform 
flows found in the high-density limit ( [Anderson fc Jackson 199*2 ). 



Typical solutions {vo{z),uq{z),Tq(z)) for steady uniform flows are given in figure ^. 
They are obtained for 9 = 20.5° by increasing the mean density i> i.e. along the thin 
solid line 9 = 20.5° in figured. As the mean density is increased, the shear close to the 
plane increases (see figure pp). The difference of temperature between the plane and 




Figure 3. (a) Density profiles vo{z), (b) velocity profiles uo{z) and (c) temperature profiles 
To{z) for steady uniform flows at a fixed angle {9 = 20.5°) when the mean density is increased. 
The inset gives the value of {h, v) for the different profiles (e = 0.6, = 0.05). 



the free surface also increases with the mean density (see figure ^). Note that with 
our boundary conditions, the temperature is always higher at the plane than at the free 
surface. Finally, the density profile vq{z) is strongly modified when varying the mean 
density (see figure ||a). For dilute flows, the thickness h is large and the density profile 
is non-inverted i.e. the density decreases with z. As the mean density increases, the 
flow becomes thinner and the free surface is better defined. At a given mean density, 
the density profile becomes inverted i.e. the maximum of density is no longer at the 
plane but in the bulk. In the {h, v) diagram of figure ^, the limit between non-inverted 
density profiles and inverted density profiles is given by the solid line. Above this line, 
the maximum of density is in the bulk. 

The shape of density profiles can be understood qualitatively by inspecting the equa- 



tions (p.5|)-(3.7). Eq(3.5) shows that the sign of the density gradient results from a 
competition between gravity, which tends to increase density close to the plane, and the 
negative temperature gradient, which tends to decrease density close to the plane. The 
temperature gradient results from the shearing at the base, which heats the material 
from below and the inelastic dissipation in the bulk, which cools the material. Therefore, 
inverted density profiles are observed when heating due to the shear at the base and the 
cooling due to dissipation in the bulk create a temperature gradient strong enough to 
overcome the gravity. 

3.3. Influence of the parameters 

The results presented so far have been obtained using e — 0.6 and cj) — 0.05. Here we 
discuss briefly the qualitative influence of the inelastic coefficient e and the boundary pa- 
rameter (j) on the solutions presented above. First, inverted density profiles are observed 
only for small values of (j). For cj) > 0.3, no inverted density profile can be obtained, 
neither by changing the control parameters {h, v) nor by modifying the coefficient of 
inelasticity e. For high values of the plane is indeed no longer a source of fluctuating 
energy (see equation 3.12| ). The density gradient close to the plane is therefore negative 
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and no inverted density profile can be obtained. For values of less than 0.3, the results 
are qualitatively the same as for (f> — 0.05. Varying changes the range of inclinations 
where steady uniform flows are observed. For example, the maximum angle is increased 
from 6 ~ 23° to 9 ^ 32° when increasing^ from 0.05 to 0.15 (e = 0.6). No qualita- 
tive change is also observed when changing the inelasticity e. Increasing the value of 
e decreases the angles and shrinks the domain of steady uniform flows. With e = 0.8 
(</) = 0.05), the range of angles is 11° < 6* < 18°. 

Steady uniform flows with inverted density profiles can therefore be obtained in the 
framework of the kinetic theory of granular flows. The question we address in this study is 
the stability of these inverted density profiles under transverse perturbations. Obviously, 
gravity is a destabilizing effect since the heavy fluid is above the light one. However, 
gravity has to overcome the stabilizing effects due to viscosity and thermal conductivity. 
In a rapid granular ffow, viscosity, thermal conductivity and density proflles are coupled 
to the flow and the prediction of the stability properties is not straightforward. 



4. Three-dimensional linear stability analysis 

We investigate the stability of steady uniform flows (vq, uq Tq) found in the previ- 



ous section using the classical normal mode analysis (Drazin & Reid 1981). The basic 
flow is perturbed by inflnitesimal disturbances, and their time evolution are studied by 
linearizing the governing equations about the basic state. The perturbations are then 
decomposed into different Fourier modes and because of the linearity of the governing 
equations, the stability of each mode can be analyzed separately. 

4.1. Normal mode analysis 

As point out by Alam & Nott (1998), the Squire theorem does not hold for rapid granular 
flows since density and temperature are not constant across the layer. Therefore, the 
first instability is not necessary two-dimensional. In this study, we are interested in 
the formation of longitudinal structures and we restrict our analysis to flows which are 
invariant along the x-direction. The flow {v, u, v, w, T) is perturbed around the basic 
flow (i^Oj Wo, To) a-iid written as: 

{u = Mo(z) + ui{y,z,t) 
v = vi{y,z,t) , T = Toiz) + Ti{y,z,t), (4.1) 

w = wi{y,z,t) 

where vi, (mi,wi,wi) and Ti are respectively the density, velocity and temperature dis- 



turbances. By substituting (4.1) into the governing equations (pTl|)-(2.3) and then lin- 
earizing about the basic state, we obtain a set of linear equations for (j^i, mi, wi, wi, Ti) 
(see Appendix A.l). Then, we seek normal mode solutions for density, velocity and 
temperature perturbations: 

(i^i,Hi,«i,u.i,ri) = 5Re [(^v{z),u{z),v{z),w{z),f{z)) e"*+*'=^] . (4.2) 

We have restricted the stability analysis to temporal stability i.e. the growth rate a is 
complex whereas the transverse wavenumber k is assumed to be real. The basic flow 
(fo, uo, To) is unstable under the transverse perturbation of wavenumber k if the real 
part of the growth rate, 5Re(o'), is positive. After some algebra, it can be shown that the 
perturbed variables satisfy a system of ordinary differential equations: 

Lo{z)^X{z) + Mo(z)-^X(z) + Nq{z)X{z) = 0, (4.3) 
dz^ az 
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where X{z) is the five-elements vector defined by X{z) — {£'{z),ii(z),v{z),w{z),T{z)), 
and Lo{z), Mo{z), No{z) are 5x5 matrices which are given in Appendix A. 2. Note that 
these matrices depend on the basic flow (i^q, uq, Tq), on the wavenumber k and on the 
growth rate a. The boundary conditions at the plane can be written in the same manner: 

Qo{z)-^X{z) + Roiz)X{z) = atz = 0, (4.4) 
dz 

where Qo{z) and Roiz) are matrices which are also given in Appendix A. 2. At infinity, the 
disturbances X(z) have to vanish. It can be shown by doing an asymptotic expansion of 
equations ( |4.3| ) that the disturbances (m, v, w, T) decrease as exp(— fcz) when z is much 
larger than the characteristic thickness of the basic flow (see Appendix B). Matching 
the boundary conditions at infinity thus leads to numerical difficulties when k becomes 
small since the computational domain varies as 1/fc. Instead of writing the boundary 
conditions at infinity, the equations are integrated up to a finite value z = Zmax and the 
disturbances (w, it, w, T) are assumed to satisfy 

^ [u{z),v{z), W{z), f{z)^ = -k {u{z), V{z),w{z), f{z)^ , at 2; = Z^nax- (4.5) 

The system of five ordinary differential equations ([4.3[) together with the eight boundary 



conditions {iA) and ( |4.5| ) constitute an eigenvalue problem. For a given basic ffow (i/Q, 
Wo, Tq) and wavenumber fc, a non-zero solution (i>, m, w, w, T) exists only for specific 
values of the growth rate a. 



4.2. Numerical method 

We have solved the above eigenvalue problem thanks to a Chebychev spectral collocation 
method (Gottlieb et al 1984). Chebychev collocation approach has shown to be well- 
adapted to the stability of boundary-layer ffows, since Chebychev polynomials resolve 
the boundary regions extremely well ( Malik 199C| ) . It is then suitable for our ffow which 
is localized close to the plane. Moreover, the use of collocation makes the derivatives 
easy to compute and simplifies the treatment of boundary conditions. 
The principle of the Chebychev spectral collocation method is to discretize the ordinary 
differential equations ( |4.3| ) by interpolating the perturbed functions {v{z), u{z), i{z), 
w{z), T{z)) on + 1 collocation points given by 



cos 



TTJ 

N' 



(z = 0,...,iV). 



(4.6) 



In order to relate the Chebychev space (<j e [—1, 1]) to the physical domain (z e [0, Zmax]), 
we use a two-parameters algebric transformation ( |Malik 1990 ) 



1 + ^ 



(4.7) 



where 6=1-1- 2a/zmax and a = zuZmax I {zmax — 2z/i). The location Zh corresponds to 
C = 0, i.e. half of the grid points are located in the region Q < z < Zh- This mapping 
allows to cluster grid points near the plane. Using the expression for the derivatives at 
the collocation points (Gottlieb et al 1984), the system of ordinary differential equations 
(4.3) together with the boundary conditions (4.4) and ( |4.5| ) reduce to a linear algebraic 
eigenvalue problem that can be written in the following form: 



AX = aBX. 



(4.8) 
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Here, X is the vector containing the 5(iV+ 1) elements of the interpolation of {iy{z), u{z), 
v{z), 'w{z), T(z)) on the {N + 1) collocation points: 



,un; Wo, wo,...,wn] Tq,...,Tn 



(4.9) 



and A and B are b{N + 1) x 5(7V+ 1) matrices computed from the matrices Lq, -^^o and 
Z/Q. A and B depend on the basic flow and on the wavenumber k. Since the boundary 
conditions (4.5) do not contain the eigenvalue ct, the matrix B is singular. It must be 
noticed that the above discretization requires ten boundary conditions while the physical 
boundary conditions ( [f.4| ) and ( |4.5| ) give only eight relations for the perturbed fields 
(four at the plane and four at infinity) . The same problem arises in the stability analysis 
of compressible hydrodynamic flows (Malik 1990). Here we have chosen the vertical 
momentum balance at the plane and the conservation of mass at infinity as the two 
extra boundary conditions. 



The generalized eigenvalue problem (4.S) is solved by using the scientific software Mat- 
Lab. The advantage of spectral methods is that the whole spectrum of eigenvalues may 
be ob tained. However, ma ny eigenvalues are spurious eigenvalues due to the discretiza- 
tion ( Mayer fc Powel 1992 ). The location of these spurious modes are very sensitive to 
the number N of collocation points while the physical modes are insensitive. This allows 
us to select the few physical modes among the large spectrum of the discretized problem. 
With a typical number of collocation points = 50, the absolute error in the physical 
eigenvalues is less than 10~^. 



5. Results 

5.1. Diagram of stability 

Figure ^ gives the diagram of stability in the plane (/i, v) for e = 0.6 and (j) — 0.05. The 
region of non-inverted density profiles is presented on the same diagram. Since we are 
essentially interested in thin granular layer flows, we have limited our stability analysis 
to /i i 30 particle diameters. The study is also limited to mean density less than 0.4. 
For higher mean density, the density profile is very sharp and the number of collocation 
points required for accurate computation increases dramatically. 

The stability diagram can be divided in three regions: a stable region (hatched region), 
an unstable region with stationary modes, an unstable region with propagating modes. 
The stable zone is delimited by the marginal curve where the real part of a vanishes. One 
observes that for h less than 7 particle diameters, the flow is always stable whatever the 
value of the mean density. The unstable region is composed of two regions: one where 
the most amplified mode is stationary i.e. 5m(CT)= 0, and one where the most amplified 
mode is propagating i.e. with a non zero phase velocity (3m((T)7^ 0). Propagating modes 
are obtained at low density while stationary instability takes place at high density. 

The transition between the stationary instability and propagating modes appears 
clearly in figure ||. We have plotted in figure ^ the growth rate 3fie(cr) and 3m(cr) as 
a function of the wavenumber fc, for the most dangerous mode at different locations in 
the stability diagram (A, B and C in figure We keep the thickness constant {h = 12.8) 
and increase the mean density from v = 0.08 to v = 0.16. One sees that the station- 
ary mode in C appears from the collapse of the two conjugated propagating modes. In 
figure ^ we also observe that at the threshold of the instability, the first unstable mode 
occurs at a finite wavenumber kc- We have systematically studied the critical wavelength 
Ac = 27r/fcc along the marginal curve. When the instability is stationary, the wavelength 
Ac is nearly constant and equal twice the thickness of the flow. By contrast, the wave- 
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Figure 4. Stability diagram for e = 0.6 and 4> ~ 0.05. Tlie flow is stable inside the hatched 
region and unstable outside. The bold solid line is the marginal curve. The dotted line separates 
a region where the most amplified mode is propagating from a region where it is stationary. The 
grey zone is the domain of non-inverted density profiles. 
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Figure 5. Growth rate 5Re((T) and 9m((T) as a function of the wavenumber k for the most 
dangerous mode. The thickness is constant {h — 12.8) and the mean density is increased: 
i> = 0.08 (A), = 0.13 (B) and v = 0.16 (C) (see the stability diagram in figure |). e = 0.6 and 
(h = 0.05. 



length Ac strongly varies along the marginal curve for the propagating instability. In 
figure ^ we have plotted Xc/h as a function of h. We observe that Ac is about 2.5 the 
thickness h of the flow as long as ft. < 15 particle diameters but starts to increase for 
large thickness. For example, h — 21 correspond to a wavelength Ac = 315, which is 15 
times the thickness of the flow. We shall discuss later this behavior as the signature of a 
change in the instability mechanism. 

5.2. Eigenf unctions 

The three-dimensional stability analysis of the basic flow gives also the five eigenfunctions 
{v{z),u{z).,v{z),w{z).,T{z)) for a given wavenumber k and growth rate a. In figure 0, 
we have plotted the real part of the eigenfunctions corresponding to the most amplified 
mode for (/i = 8, P = 0.23). This mode is stationary. We observe that the perturbation 
of density i'(z) vanishes for z > h whereas the velocity and temperature perturbations 
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h 

Figure 6. \c/h (•) as a function of h , where Ac is the wavelength selected by the instability 
along the marginal curve for the propagating instability. The white circles (o) is the same curve 
obtained when the coUisional dissipation 7 is artificially set to zero, (e = 0.6, = 0.05). 

extend further and decay exponentially (see appendix B). From these eigenfunctions, we 
can simply recover the perturbed field {vi, ui,vi, wi,Ti){y, z, t). For example, viiy, t) 
is obtained from the eigenfunction v{z) by the relation vi(y,z,t) = 3fie[i>(z) exp(crt + 
iky)] (see equation 4.2). In the following the perturbed field are plotted for t = 0. In 
figure ^, we present the perturbations of the flow corresponding to the most amplified 
stationary mode presented in figure ^ (ft, = 8, P = 0.23). The pattern is typical of that 
observed in the stationary unstable region. Figure |^(a) presents the transverse velocity 
field {vi (y, z, 0), wi{y, z, 0)). We can see that the motion in the transverse plane consists 
in counter-rotating vortices with a pair of vortices per wavelength A. This transverse 
flow is strongly correlated with the perturbation of the longitudinal velocity. Figure ^(b) 
shows the contours of constant value for the perturbed longitudinal velocity mi(?/, z,0). 
The material going towards the plane is flowing faster in the slope direction than the 
material rising up to the free surface. Therefore, the 3D structure of the perturbed flow 
consists in longitudinal vortices with transverse variation of the longitudinal velocity. 

The corresponding perturbations of density vi{y, z, 0) and temperature Ti{y, z, 0) are 
given in figure ||(c) and (d). The perturbed density results from two effects: the advection 
of the basic flow by the transverse perturbed flow and the compressibility. In order to 
compare more easily the density with the transverse flow, we have plotted in figure ^ the 
depth averaged density /x(y) = /q^"°^ '^liv, z)dz as a function of the transverse direction 
y (solid line). The vertical velocity is given on the same plot (dotted line). We observe 
that the average density is higher where the flow is going downwards and smaller when 
the flow is going upwards. The perturbation of density can also be interpreted in term 
of free surface perturbation. If we define the free surface as being the surface where the 
total density (1^0(2) + i^i(j/, z, 0)) is constant and equal to 1% of the maximum value of 
Po{z), we obtain the bold sohd line in figure ||(a). One can see that the downward part 
of the flow corresponds to a through while the upward part of the flow corresponds to a 
crest. 

Most of the features of the stationary modes are recovered for the propagating modes 
except that the whole pattern is now drifting in the transverse direction, due to the non- 
zero phase velocity (see figure |l^). Qualitative differences exist since the propagating 
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Figure 7. Real part of the eigenfunctions {C'{z),u{z),v{z),w{z),T{z)) corresponding to the 
stationary mode {h = 8, u = 0.23, k = 0.4). (e = 0.6, = 0.05). 

modes lead to a phase shift between the eigenfunctions. The vortices are asymetric as 
shown in figure |lC|(a) and the transverse variations of density and longitudinal velocity 
are no longer in phase but slightly shifted. 

5.3. Influence of the parameters 

The results of the stability analysis presented in the previous section have been performed 
using a given set of parameters e = 0.6 and </> — 0.05. Changing these parameters modifies 
the solutions for steady uniform flows and therefore modifies the results of the stability 
analysis. However, the main results of the previous section are not qualitatively changed 
in the range of parameters (e, (f)) where inverted density profiles are observed. For 
instance, we have perform the stability analysis using e = 0.8 and (j) = 0.12 and recover 
the instability and the formation of longitudinal vortices. The domain of propagating 
instability increases with higher values of (f> and the vortices are more tilted from the 
vertical for the propagating unstable modes (see figure pT|) . 

However, an important difference with the case (e = 0.6, — 0.05) is that the curve of 
marginal stability may overlap the domain of non- inverted density profiles. This means 
that for some values of 4> and e, it is possible to find basic flows with non-inverted density 
profiles which are unstable. These non-inverted unstable flows are observed in a narrow 
range of the stability diagram for large thickness {h > 30) and small density {D < 0.1). 
We will see in the next section that inelasticity is responsible of the instability of such 
non-inverted density profiles. 

6. Discussion 

6.1. Instability mechanism and role of the inelasticity 

In figure |[ the domain of non-inverted density profiles is shown together with the domain 
of instability. The correlation between the two regions suggests that inversion of the den- 
sity profile plays an important role in the instability mechanism. We have checked the 
role of the gravity by artificially modifying g in the linearized equations. In a wide range 
of thickness and mean density, increasing gravity increases the growth rate whereas de- 
creasing gravity stabilizes the flow. In that case, the instability comes from the inversion 
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Figure 8. Perturbed fields corresponding to the stationary mode {h = u = 0.23, k = 0.4); 
(a) transverse velocity field, the bold solid line gives the free surface deformation; (b), (c), (d) 
contours of constant values for longitudinal velocity, density and temperature. The grey level is 
white for the most positive value and black for the most negative value, (e = 0.6, (j) = 0.05). 



of the density profile, which is due to the sclf-hcating at the plane ( "Rayleigh-Benard" 
type of instabihty mechanism). In fluid mechanics, the Rayleigh-Benard instability is 
controlled at the threshold by a single dimensionlcss parameter, the Rayleigh number 
Ra = gpAph^ /r]K. In a granular flow, it is difhcult to define such a non-local control 
parameter since the Boussinesq approximation is far from being satisfied. Indeed, flow 
quantities (density, temperature) as well as transport coefficients (viscosities, conductiv- 
ity) strongly vary inside the flow. 

As pointed out in the previous section, some flows with non-inverted density profiles are 
unstable (e.g. with e = 0.8 and (j> = 0.12). This means that gravity is not the only desta- 
bilizing effect. Another well-known source of instability in granular flows is inelasticity. 
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Figure 9. Perturbed depth averaged density fi{y) = JJ'""'" vi{y,z)Az (solid line) as a function 
of y for the stationary unstable mode (/i = 8, = 0.23, k = 0.4). The dotted line presents 
the transverse variations of the perturbed vertical velocity wi (arbitrary scale). The perturbed 
averaged density is higher where the flow is going downwards {w\ < 0) and smaller when the 
flow is going upwards {wi > 0). (e = 0.6 and (p — 0.05). 



Studies on two-dimensional shear flows have shown that the dissipation due to inelastic 



collisions can lead to the formation of clusters (Savage 1992 Chi-Hwa Wang et al 1997 



Alam & Nott 1998; Nott et al 1999). In order to better understand the role of dissi- 
pation in our problem, we have performed the stability analysis by setting to zero the 
collisional rate of energy dissipation 7 in the linearized energy equation. In all the other 
terms, inelasticity is kept in order to study the same basic flow. Figure |l^ presents on 
the same diagram the curve of marginal stability in the case 7 7^ and 7 = (e = 0.6, 
(j) = 0.05). We observe that the flow remains unstable in a wide range of parameters 
with 7 = 0. This proves that inelasticity is not necessary to get an instability. For 
thin and dense flows, inelasticity stabilizes the flow whereas for dilute flows inelasticity 
slightly lowers the threshold. The interesting point is that non-inverted density profiles 
are always stable without inelasticity, whatever the choice of parameters e and cf). The 
dissipation also strongly influences the wavelength selection Ac at the threshold. In flgure 
^ we have plotted Xc/h as a function of h both for 7 7^ (black circles) and 7 = (white 
circles). When the dissipation is zero 7 = 0, the wavelength at the threshold scales with 
the thickness of the flow over the whole range of thickness (Ac ^ 2h). This implies that 
the increase of Ac for large h observed in the real system (7 ^ 0) comes from inelasticity. 

This analysis suggests that both gravity and inelasticity contribute to the instability 
depending on the parameters. For thin and dense flows, gravity is the dominant desta- 
bilizing effect giving rise to a Rayleigh-Benard kind of instability, whereas for large and 
dilute flows inelasticity is the principal ingredient of instability, giving rise to a clustering- 
like instability. However, it is not possible to separate the two effects since gravity and 
inelasticity are associated to the same unstable mode. 

6.2. Comparison with experimental observations 

The present study captures the principal features of the longitudinal vortices instability. 
First, our analysis has shown that steady uniform flows down inclined planes may be 
unstable under transverse perturbations. In the range of parameters of the experiment 
{h = 10 — 12 grain diameters, P = 0.2 — 0.3), the physical origin of this instability is 
the inversion of the density profile, which is induced by the self-heating at the plane. 
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Figure 10. Perturbed fields corresponding to the propagating mode {h = 10, = 0.13, k — 0.2); 
(a) transverse velocity field; (b), (c), (d) contours of constant values for longitudinal velocity, 
density and temperature, (e = 0.6, (j) = 0.05). 

Secondly, the unstable flow consists in longitudinal vortices leading to surface deforma- 
tion, transverse variations of longitudinal velocity and density in agreement with the 
experimental observations. As observed in the experiment, the longitudinal velocity is 
larger in the troughs (i.e. where the flow is going downwards) than in the crests (i.e. 
where the flow is going upwards). The same variation occurs with the density, troughs 
are dense and crests are dilute. Finally, the instability selects a wavelength which is 
always 2 — 3 times the average thickness of the flow above the threshold, as observed in 
the experiment. 

It is more difficult to conclude about the relevance of this analysis in order to explain 
the phase velocity observed in the experiment. We have seen that depending on the 
control parameters {h, 9), the growth rate of the instability may be complex i.e. ai ^ 0. 
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Figure 11. Typical transverse velocity field with e — 0.8 and (/> = 0.12 for the propagating 
instability {h = 17, v = 0.25, k = 0.3). 




Figure 12. Role of inelasticity on the stability diagram. The solid line is the marginal curve 
with 7 5^ and the dotted line is the marginal curve with 7 = (e = 0.6, = 0.05). The grey 
zone is the domain of non-inverted density profiles. 



This imaginary part of the growth rate is associated to a phase velocity — ai/k, where 
k is the wavenumber. A typical order of magnitude for is (0.1 — 0.5)-y/gd which is 
a few percents of the chute velocity. This order of magnitude is compatible with the 
drift velocity observed in the experiment. However, the theory predicts a transition be- 
tween propagating and stationary instability, which is not observed in the experiment. 
Moreover, in the range of parameters of the experiment, the theory predicts a stationary 
instability. The phase velocity observed experimentally is more likely related to the very 
weak inclination of the pattern in the experiment. This inclination is not taken into 
account in this study since we have restricted our analysis to pure transverse pertur- 
bations {kx =0, ky ^0). It would be interesting to investigate the stability of steady 
uniform flows when a small perturbation in the longitudinal direction {kx << 1) is added 
to the ky perturbation. This should give rise to a phase velocity in the whole range of 
parameters. 

The agreement between theory and experiment is only qualitative. The most impor- 
tant discrepancy between theory and experiment lies in the range of angles where the 
instability is predicted. With e = 0.6 and (j) = 0.05, the predicted angles are ~ 20° 
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while in the experiment the minimum angle in order to observe the instability is 9 — 38° 
with sand 0.25 mm diameter. The same order of magnitude is obtained in the exper- 
iment with monodisperse glass beads. In the theory, it is not possible to reach such 
inclination by varying the parameters (e, 0) except for non-physical values of inelasticity 
(e ~ 0). Actually, this discrepancy between theory and experiment is not surprising since 
the standard kinetic theory is known to be limited. The kinetic theory has been first 
developed for quasi-elastic particles (e ~ 1) and rather dilute flows. Even under these 
conditions, experimental studies of 2D coUisiona l granular flows d own inclined planes 
have shown that the theory is not quantitative ( Azanza et al 199E ). One of the main 
problem of the kinetic theory is the determination of the boundary conditions at the 
rough plane. Experiments in dicate that, even in the coUisional regime, the material 
is structured near the plane ( Azanza et al 1999 ). This local organization is not taken 
into account in the kinetic theory. Since boundary conditions control the production 
of granular temperature, it is not surprising that the theory fails in quantitatively de- 
scribe steady uniform flows. Another difficulty of the kinetic theory lies in the treatment 
of strongly inelastic particules. For high inelasticity, additional terms should be intro- 
duced in the kinetic theory in order to take into account the lack of separation between 



the microscopic and macroscopic scales inherent to inelasticity (Tan & Goldhirsch 1999 



; Sela & Goldhirsch 1998). However, such additional terms complicate considerably the 
equations and are still the subject of active researches ( poldhirsch 1999 ). 



7. Conclusion 

In this paper, we have performed a three-dimensional stability analysis of rapid gran- 
ular flows in the framework of the kinetic theory of granular gases. We have shown that 
steady uniform flows down inclined planes can be unstable under transverse perturba- 
tions. The structure of the unstable flow consists in longitudinal vortices with transverse 
variations of free surface, chute velocity and density in agreement with the experimental 



observations (Forterre & Pouliquen 2001). Actually, the agreement is only qualitative. 
This is not surprising since the experiments take place in a semi-dilute regime and use 
rather inelastic particles, which is beyond the domain of applicability of the simple ki- 
netic theory used in the analysis. However, our study shows that the kinetic theory is 
a relevant framework for the description of rapid granular flows. The kinetic theory has 
revealed a new instability mechanism specific to granular material: inelastic collisions 
trigger a self-induced convection yielding longitudinal vortices in chute flows. 

In classical fluid mechanics, Rayleigh-Benard convection is the paradigm for pattern 
forming instabilities and the study of the transition to turbulence. Therefore, an impor- 
tant question is whether the granular convection observed in our experiment represents 
the starting point of a similar scenario towards more complex structures. Actually, 
we have observed in the experiment a new pattern when the plane is strongly inclined 
{9 > 50° with sand 0.25 mm in mean diameter). Instead of longitudinal streaks, a regular 
square pattern looking as "fish scales" develops on the free surface, as shown in figure 
Note that under such inclinations, a fast camera is necessary to capture this structure. 
The appearance of this new pattern raises several issues. A first possibility is that for this 
range of parameters the square pattern represents the most unstable mode for the pri- 



mary instability, as for example observed in inclined layer convection (Daniels et al 2000). 
This could be investigated by generalizing the present stability analysis to longitudinal 
modes i.e. and ky 0. A second possibility is that the "scales" result from a 

secondary instability of the longitudinal vortices, as a consequence of nonlinear effects. 
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Figure 13. Top view of the free surface of the flow showing the formation of "scales" when the 
plane is strongly inclined (sand 0.25 mm mean diameter, = 52°, hg = 15 mm). The picture is 
taken with a short shutter time of 1/lOOOOs. 



Such an evolution is well-documented for classical fluid flows (Godreche & Manneville 
1998) but still remains an open issue for granular flows. 
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Appendix A. 

In this Appendix we detail the main steps of the linear stability analysis presented in 
Section 4.1. 

A.l. linearization of the governing equations 



The linearisation of the mass equation (2.1), the momentum equation (|2.2|) and the 



energy equation ( |2.3| ) around the basic flow (i^q, uq, Tq) gives: 

di^i f dvi dwi \ duo , . 

+ ^^^^ 
dui , duo \ , „ , dY^xy^ OY^xzi 
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where Si is the stress tensor disturbance, qi is the heat flux disturbance and 71 is the 
energy dissipation disturbance given by: 
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Here the pressure disturbance Pi, the viscosity disturbance iji, the conductivity distur- 
bance Ki and the energy dissipated disturbance 71 are hnearized functions of i>i and Ti 
which are given in table 2. 

must be considered together with the boundary con- 



The system of equations ( A 1 )- 
ditions for the disturbances {vi{z), wi{z), Ti{z)). At the plane, the boundary 

conditions are written as: 
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at z = 0, 
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at z = 0, 
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at z = 0, 



(A 6) 
(A 7) 

(A 8) 
(A 9) 



where rji* and Ki* are given in table 2. The first condition (A 6) simply expresses that 
t he p lane i s rigi d. The last three conditions (A7)-(A9) are the boundary conditions 
( ^T0| ) and jlnl ), which are linearized around the basic flow (i^q, uq, Tq). At infinity, we 
assume the disturbances to vanish: 



vi,ui,vi,wi,Ti 0; 



when , 



(A 10) 



A. 2. Matrices of the eigenvalue problem 



The non-zero elements of the coefficient matrices in equations (4.3) and (4.4) are given 
below. The definition of the functions of the basic fiow (ao,. . . j^o) ai'e given in table 2. 
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-Pi = fi' i'^o)Tai^i + fi{i^a)Ti = a^i'i + fooTi 

VI = /2'(i'o)ro5z/i + if2{uo)To-^Ti = cofi + doTi 

Ki = fs'MTo^i^i + i/3(i'o)To-iTi = eoi^i + hoTi 

71 = fs'MToiui + |/5(i'o)Toiri = loui+moTi 

Vi* = <A/6'(i'o)Toii/i + ^(pfeMTo'^Ti = hqi^i +PoTi 

71* = (1 - e')fr'{uQ)Tohi + ^(1 - e')/7(i/o)To'ri = qoi^i + roTi 

Table 2. Pressure disturbance Pi, viscosity disturbance iji, conductivity disturbance Ki and 

energy dissipated disturbance 71 
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Appendix B. 

Here we give the asymptotic behaviour of the disturbances {v, it, v, w, T) when z is 
much larger than the characteristic thickness of the basic flow. The linearized equations 
(4.3) which govern the disturbances can be written as 



-vp ikv H — — I — w 



, .duo J 

1^0 [ (TU+ 



VqGV 



Bl) 
(B2) 

(B3) 



VpGW = 
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dTo 
dz 
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(B4) 



k^KoT 



(B5) 



In order to obtain the asymptotic expression of these equations, one has to know the 
asymptotic behaviour of the basic flow. Equation ( |3.5| ) shows that the zero-order density 
Uq decays exponentially to zero as 



exp 



To 



when z is large. Using equations (3.5)-(3.7), the asymptotic behaviour of the basic flow 
is given by 
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Considering (B 6)-( B 11 ) into the linearized equations (B2)-(B5) gives an asymptotic 
behaviours for z which is much larger than the characteristic thickness of the basic flow. 
From equation (BJ), one shows that the perturbed density decays on the same length 
scale as the basic flow: 

vr^vo. (B12) 
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Therefore, the asymptotic hmit of equations (B 1)-(B 5) is given by 



These equations are easily solved and gi 




(B13) 
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